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TECHNIQUES  FOR  DIRECTIONAL  DATA 
by 

M.A.  Stephens 

McGill  University,  Montreal  and  Stanford  University 

1.  INTRODUCTION 

1.1  In  recent  years  techniques  have  been  developed  for  dealing 
with  statistical  data  where  the  observations  are  directions,  and  where 
the  directions  are  assumed  to  be  more  or  less  concentrated  around  a 
single  mode*  In  three  dimensions,  the  distribution  used  to  describe 
such  directional  data  is  Fisher's  (1953)  distribution,  and  in  two 
dimensions  it  is  the  von  Mises  distribution.  In  this  paper  we  extend 
the  techniques  for  these  distributions  to  deal  with  axial  data,  i.e. 
data  consisting  of  vectors  whose  direction  can  be  in  either  sense,  and 
also  for  use  with  directed  data  from  populations  with  two  modes,  in 
opposite  directions.  The  techniques  make  use  of  tables  already  prepared 
for  the  Fisher  and  von  Mises  distributions-  Examples  of  directed  data 
are,  in  three  dimensions,  directions  of  magnetization  of  rocks,  or,  in 
two  dimensions,  directions  of  bird  flights  or  of  prevailing  winds; 
examples  of  axial  data  are  normals  to  planes  of  cleavage  of  locks,  or 
inclinations  of  she  long  axis  of  pebbles  in  till  deposits. 

The  procedures  will  be  given  first  for  three  dimensions,  since 
later  on  it  is  easy  to  adapt  them  for  two  dimensions;  the  rest  of  the 
introduction  deals  with  notation  to  be  used.  The  Fisher  distribution  is 
described  in  section  2;  its  bimodal  extension,  assuming  axial  data  or 
directed  data  with  equal  modal  strengths,  is  given  in  section  3,  with 
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examples}  tests  of  hypotheses  for  three-dimensional  data  are  In  section 
U.  For  two  dimensions,  the  von  Mises  distribution  and  Its  blmodal 
extension  are  treated  in  section  5 • 


The  adaptation  of  all  the  procedures  for  directed  data  but  with 
unequal  modes  in  given  in  section  6.  Section  7  demonstrates  the  impor¬ 
tance  of  knowing  the  type  of  data,  and  deals  briefly  with  related  topics. 
Examples  are  included  throughout  the  paper. 


1.2  Notation.  Observations  denoting  directions  sure  recorded  by  unit 
vectors}  in  three  dimensions  a  typical  vector  is  OP,  starting  at  the 
center  0  of  a  sphere  of  radius  1  and  ending  at  P,  a  point  on  the 
surface  of  the  sphere.  Axial  data  could  be  best  recorded  by  drawing 
the  entire  diameter  POQ  say,  though  in  three  dimensions  this  is  then 
difficult,  in  practice,  to  show  on  the  usual  diagrams.  Techniques  for 
axial  data  must  not  depend  on  whether  gP  or  OJ  is  used  to  represent 
an  observation}  for  all  data,  thus  the  vector  used  is  called 
OJ5  .  Let  P  be  located  by  the  usual  spherical  polar  coordinates  9, 

41;  we  shall  regard  these  also  as  coordinates  of  the  line  OP.  For 
simplicity,  let  6=0°  be  thought  of  as  pointing  to  the  "north  pole" 
of  the  sphere,  so  that  9  =  90°  is  the  equator  and  6  =  l80°  is  the 
"south  pole". 


1-3  Other  Notation.  In  the  northern  hemisphere,  6  is  then  the 
co latitude  of  P,  and  9  =  90°-  X,  where  X  is  northern  latitude} 
in  the  southern  hemisphere,  9  -  90° +  X,  where  X  is  southern  latitude; 
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<t>  is  the  longitude  measured  from  0  to  360°  eastwards  from  <t>  «  0°. 

<t  is  called  also  the  orientation}  P  is  sometimes  measured  by  orienta¬ 
tion  4>  and  dip  angle  5  below  the  equator}  6  corresponds  to  southern 
latitude  X,  and  0  =  90°'*'X..  In  practical  work,  0  (or  6)  and  #  are 
often  given  in  degrees,  as  has  so  far  been  assumed*  in  theoretical  dis¬ 
cussion,  we  shall  assume  0,  ®  are  in  radians.  This  will  not  affect 
the  practical  techniques,  which  use  the  components  of  the  given  vectors. 
For  these,  we  must  introduce  a  set  of  rectangular  coordinates}  a  natural 
set  has  the  z-axis  along  0=0,  the  x-axis  along  0  =  rt/2,  4>  =  0, 
and  the  y-axis  along  9  =  n/2,  $  =  n/2.  The  components  of  a  vector 
with  coordinates  0,  ®  are  then 

x  =  sin  0  cos  $  ,  y  =  sin  0  sin  4>  ,  z  =  cos  9  . 

For  a  given  sample  of  size  N,  let  £1^  (i=l,2, . . .,N)  be  the 
i-th  unit  vector,  and  let  x^  y^  z  be  its  components)  define 
X  =  Ex^,  Y  =  Ey^,  Z  =  Ez^  .  These  are  the  components  of  the  vector  sum 
or  resultant  g  of  the  ft^}  if  g  has  length  R,  then 
R2  =  X2  +  Y2  +  Z2.  When  dealing  with  several  samples,  the  subscript 
r  will  give  the  value  for  the  r-th  sample,  e.g.  Nr  is  the  sample 
size  of  the  r-th  sample,  with  resultant  g^  N  will  denote  the  total 
sample  size  +•••+  Ng,  where  s  is  the  number  of  samples. 

2.  TEE  FISHER  DISTRIBlfTION  FOR  DIRECTED  VECTORS. 

2.1  Suppose  the  vectors  Oj^  (1=1,2, . . .,N)  represent  directed 
data,  with  an  arrowhead  at  P  . 


3 


The  Fisher  distribution  describes  a  probability  density  at  the 
point  9 <t>o  which  is  proportional  to  exp(*c  cos  0q):  more  precisely, 
the  density  of  6  is 

(l)  f(9)  =  exp(<  cob  8)  ,  0  <9  <  n  , 

and  4>  is  independently  uniformly  distributed  between  0  and  2rt. 

The  distribution  is  symmetrical  around  &  (along  9=0,  pointing  to 
the  north  pole),  and  with  i  single  mode  at  A  ;  <  is  a  parameter 

(<  >  0),  which  describes  the  concentration  of  the  distribution.  When 
k  is  large  the  distribution  is  highly  concentrated  around  A,  end 
when  k  «  0  the  vectors  (i.e.  the  points  P)  are  uniformly  distributed 
over  the  surface  of  the  sphere-  This  distribution  was  introduced  by 
Fisher  (l953)  to  describe  vectors  denoting  remanent  magnetization  of 
rocks.  Statistical  procedures  were  given  by  Fisher,  and  by  Watson  (1956), 
Watson  and  Williams  (1956),  and  Watson  and  Irving  ( 1956);  this  work  has 
been  developed,  and  the  necessary  v.ables  produced,  Dy  Stephens  (1962b, 
1967,  1969a),.  The  techniques  have  begun  to  be  used  in  applied  work, 
particularly  in  a  geological  context;  see,  for  example,  Andrews  and 
Shimizu  <(1966),  and,  for  a  wider  discussion,  with  a  long  list  of 
ref erenees,  Watson  (1968).  It  will  be  convenient  now  to  summarize  these 
procedures;  in  the  next  section  they  will  be  adapted  for  use  with 
axial  or  bimodal  three-dimensional  data. 

2 .2  Estimation  of  modal  vector  and  concentration  parameter  for 
the  Fisher  distribution.  In  the  distribution  (l)  above,  the  modal 
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direction  A  was  assumed  along  3  =  0>  in  practice  it  will  usually  not 
be  known*  and*  with  <*  must  be  estimated*  The  maximum  likelihood 
equations  for  these  estimators  are  based  on  the  statistics  g  and  X* 
from  a  sample  of  N  unit  vectors*  already  described  in  section  1* 

(a)  The  estimate  of  A  is  the  direction  of  R,  i*e**  a  unit  vector 
a  estimating  A  is  a  =  B/B  . 

(b)  To  estimate  K*  solve 

(2)  coth  K~ljK-B/lii 

if  A  is  known*  replace  R  by  X  in  this  equation.  A  table  for 
solving  (2)  is  in  Stephens  (1967)* 

3*  BIMODAL  EXTENSION  OF  THE  FISHER  DISTRIBUTION 

3*1  The  natural  extension  of  the  Fisher  distribution  to  cover 
bimodal  data  is  obtained  by  superimposing  two  Fisher  distributions  with 
opposite  modal  vectors.  If  the  line  along  9  =  0  and  9  -  n*  now 
called  the  modal  axis  A*  represents  the  direction  of  the  two  modes, 
the  density  of  9  is  then 

(3)  f(9)  =  2~slriwf  exp(*  cos  9)+U-a)exp('K  cos  9))  , 

O  <  9  <  n 

and  0  has  a  uniform  distribution  between  0  and  2n  as  before. 

The  relative  strength  of  the  two  modes  is  measured  by  the  parameter 
a*  which  lies  between  0  and  !•  When  a  =  0.5*  (3)  becomes 


2 'sin h£^  (c08h  (*  008  S))>  °<9  <«  * 


(4)  f(9)  - 

Distribution  (3)  will  be  used  to  describe  directed  data  with  unequal 
modes)  distribution  (4)  will  be  used  in  analysing  axial  data,  if  both 
ends  are  recorded.  In  general,  the  modal  axis  will  not  be  known,  and 
must  be  estimated.  When  A  1b  known,  pointing,  say,  to  the  North  pole, 
we  could  choose  to  record  axial  data  by  one  point  only,  in  the  northern 
hemisphere;  the  density  of  0  is  then 

(5)  f(e)  *  (cosh  cos  9))>  0  <  6  <  */2  . 

In  practice,  a  given  experimenter,  say  collecting  directions  of 
magnetization  of  rocks,  will  probably  show  a  natural  preference  for 
recording  the  data  in  one  sense  more  often  them  the  other,  bo  that  a 
given  sample  of  axial  data,  as  recorded,  may  look  like  directed  data 
with  two  unequal  modes  or  even  only  one  mode)  this  can  be  very  deceptive. 
It  is  therefore  important  to  know  what  type  of  data  is  in  a  given  sample 
in  order  to  decide  the  analysis  to  be  undertaken)  one  should  not  rely 
only  on  the  appearance  of  the  sample.  We  illustrate  this  point  in 
section  -7  * 

We  now  discuss  two  estimation  problems,  for  distribution  (4))  how 
to  estimate  the  direction  of  the  modal  axis  A,  when  this  is  not  known 
(so  that  it  does  not  lie  along  9=0))  and  how  to  estimate  The 

techniques  to  be  given  will  not  depend  on  which  end  of  a  diameter  is 
chosen  to  represent  an  observation  with  no  preferred  sense)  thus  for 
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this  type  of  data  one  follows  the  procedures  using  the  observations 
exactly  as  given.  This  is  also  the  case  for  data  representing  directed 
vectors  from  a  population  with  two  opposite  modes  of  equal  amplitude. 
Modes  of  unequal  amplitude  (a  ^  0.5)  are  discussed  in  section  6. 

3  *2  Estimation  of  the  Modal  Axis  A.  Suppose  a  plane  M  is  chosen 
through  0,  characterized  by  its  normal  ;  one  end  only  of  each 
recorded  axial  diameter  is  then  chosen  to  give  a  directed  vector,  such 
that  all  the  directed  vectors  lie  on  one  side  of  M.  The  resultant  R 

rv 

is  then  calculated,  its  value  clearly  depends  on  the  choice  of  M,  i-e. 
of  g.  When  jq  is  along  the  modal  diameter  A,  so  that  the  plane  M 
is  at  right  angles  to  A,  the  expected  length  R  of  R  is  a  maximum; 
and  when  A  lies  in  plane  M,  so  that  n  is  at  right  angles  to  A, 
the  expected  length  R  will  be  a  minimum.  This  suggests  that  to 
estimate  A,  we  must  find  the  plane  M  which  gives  a  maximum  R, 
and  this  is  done  iteratively  by  the  following  method,  which  applies  to 
both  axial  and  directed  data. 

(a)  Directed  data  is  recorded  by  the  end  P.  of  the  vector  OP  ; 

i  i’ 

for  an  axial  observation  choose  either  end  of  the  diameter  to  be  initially 

‘  Ir-  the  steps  which  follow,  the  direction  of  a  sample  vector  will 

sometimes  be  reversed,  and  P^  always  refers  to  the  end  of  the  vector 

wnich  is  currently  used  to  assign  it  a  direction  OP  . 

i 

(b)  If  a  good  estimate  V,  length  V,  of  the  modal  vector  exists, 
let  the  components  of  the  unit  vector  v  =  v/v  be  1,  m,  n. 
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(c)  Suppose  the  sample  vector  OP^  has  components  x^,  y^ 

for  each  sample  vector,  calculate  a  =  v  •  OP^,  i.e.  a  =  Ix^+my^+nz^) 
If  a  is  negative,  change  the  signs  of  x^,  y^^  and  z^.  This  reverses 
the  original  direction  of  OP^  and  ensures  that  it  now  makes  an  angle 
less  than  it/  2  with  V. 

(d)  For  the  final  set  of  sample  vectors,  calculate  the  resultant 
(vector  sum)  R,  length  R. 

(e)  Let  r  =  g/R  be  the  new  v,  components  1,  m,  n,  and  repeat 
from  step  (c).  When  two  successive  values  of  r  are  identical,  stop 
the  procedure-,  let  the  final  unit  vector  r  be  called  rQ>  the  line 
along  which  rQ  lies  is  the  estimate  of  the  nodal  axis  A. 

(f)  If  no  good  estimate  V  exists  in  (b)  above,  start  as  follows. 
Lot  (l,m,n)  be  (l,0,C)  and  proceed  with  steps  (c)  and  (d).  Repeat 
with  (i,m,n)  =  (0,1,0)  and  (l,m,n)  =  (0,0,1).  Of  the  three  values 
of  R  so  obtained,  choose  that  with  the  largest  length  R  as  the 
initial  V  ,  take  v  =  Fj/R  and  continue  from  step  (b)  above. 
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3 -3  Estimation  of  K »  1-  Modal  axis  known.  If  the  direction  of  the 

modal  axis  A  is  known,  and  lies  along  9  -  0;  it  is  easy  to  derive 
the  maximum  likelihood  estimating  equation  for  K •  It  is 

1  i  N 

(6)  cosh  K - -  —  7  cos  9.  tanh  (<  cos  9,)  , 

<  M  ^  i  i 

where  cos  9.  is  the  angle  between  OP.  and  the  modal  axis  A,  chosen 
1  +**  1 

to  point  In  either  direction. 

2.  Modal  axis  not  known.  When  A  is  not  known,  ve  measure  9 
from  the  estimate  of  A,  along  calculated  in  the  previous  section. 
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This  is  analogous  to  replacing  X  by  R  when  using  equation  (2)  for 
the  Fisher  distribution. 

For  each  sample  vector,  only  cos  0^  is  actually  required,  and  if 

the  components  of  r  are  1  ,  m  ,  n  >  and  the  final  components  of 
~o  o  o  o 

OP^  are  x^  y^  Zy  cos  9^  is  given  by  cos  9^  =  1^  +  m^  +  n^. 
Every  value  of  cos  0^  will  be  positive.  Equation  (6)  is  then  solved 
iteratively  by  the  following  steps.  Let  R°  be  the  final  resultant, 
length  R°,  of  the  sample  vectors,  i.e.  =  R°/R°,  and  let  M(/c) 
be  the  right  hand  side  of  (6)  for  any  K. 

(a)  The  quantity  (I  cos  0^/N  is  R°/N,  say  Yq;  solve 

cosh  K  -  1 Jk  =  Yq  to  give  an  initial  estimate  for  x.  This  may 

be  done  using  e.g.  Table  3  in  Stephens  (1967);  if  Yq  >  8,  is 

i/U-io). 

(b)  Solve  cosh  k  -  1 Jk  =  call  the  solution  Kg. 

(c)  Solve  cosh  k  -  1 Jk  =  M(Kg);  call  the  solution  k^>  etc., 

and  repeat  this  procedure;  the  sequence  for  K  converges,  and  the  limit 

A 

is  Kf  the  estimate  of  <. 

The  above  procedure  is  proved  convergent  as  follows.  First  suppose 
K*  is  the  solution  of  (6).  Since  tanh  x  <  1,  mCk^)  <  YQ  ;  and  since 
cosh  k  -  l/ k  is  monotonic  in  k,  K2<  Kl'  tanh  x  is  monotonic 

in  x,  so  M(<2)  <  so  by  repetition  of  the  argument,  <  Kg. 

Similarly  K^  <  Ky  etc.  Also  Yq  >  M(k*),  so  k ^  >  K*>  then 
M(<1)  >  m(k*),  and  so  Kg  >  K*  :  similarly  k  >  K*  for  all  i. 

Thus  the  sequence  of  solutions  K ^  is  decreasing,  bounded  below  by 
K*,  and  so  converges.  But  cosh  <  -l/*n  -  M(<n  ^),  and  if  we  take  the 
limit,  as  n  -*  »,  on  both  sides,  we  have  the  limiting  value  of  <  =  k* 
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As  will  be  seen  in  the  examples  which  now  follow,  the  technique 
converges  very  rapidly  for  k -values  of  5  or  more. 

3.4  Examples •  The  data  are  from  measurements  of  inclination  of 
till  deposits,  kindly  made  available  by  Dr.  C-  King  of  the  Department 
of  Geography,  University  of  Nottingham,  England*  there  are  4  samples  each 
of  2  observations,  measured  to  the  nearest  5  degrees.  The  effect  of  the 
precision  of  meas'.irement  is  not  considered  in  this  paper.  Table  1  gives 
the  data  and  the  steps  in  the  estimation  procedures,  for  Sample  1.  The 
table  is  divided  into  several  parts: 

(a)  The  25  values  of  0,  <t>  are  listed  first,  in  degrees.  The 
degree  symbol  is  omitted. 

(b)  Estimation  of  R°.  Unit  vectors  along  the  x,  y,  and  z  axes 
are  taken  as  starting  values  of  v,  assuming  no  initial  approximation 
known  for  the  modal  axis*  the  resultant,  length  is  given,  together  with 
its  coordinates  0,  <t>,  and  in  the  last  column  are  listed  those  vectors 
which  must  be  reversed  from  the  original  given  direction  to  lie  at  an 
angle  less  than  90  with  the  current  v. 

(c)  The  largest  R  obtained  from  (b)  is  20.20,  with  vector 
No.  1  reversed*  this  R,  reduced  to  unit  length,  becomes  the  next  v 
with  direction  cosines  1  =  -0.329,  m  =  0.9211,  r.  =  -0.208.  Now  vector 
No.  1  is  returned  no  its  original  sense  (thus  all  the  vectors,  as  given, 
are  within  90  of  the  current  v)  and  the  new  R  =  20.54*  or.  using  this 
to  make  the  new  v,  we  get  no  charge  in  R,  so  that  this  is  the  final 
resultant,  R°j  its  direction  cosines  are  -  0.4l8,  0.889,  -  0.188, 

and  its  coordinates  are  0  =  100.8,  $  =  115 .2.  Since  it  is  the  R°  for 
sample  1,  it  is  designated  R°  . 
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(d)  Estimation  of  <■  The  first  estimate  of  k>  derived  from 
R°  =  20.54,  is  =  5*6l.  The  successive  estimates  converge  rapidly 
to  k  =  5 *56. 

Table  2  gives  the  final  results  for  Samples  2,  3  and  4,  for  Samples  1 

and  2  taken  together  as  one  sample  of  50  observations*  and  for  Samples 

3  and  4  taken  together.  X,  Y ,  Z,  are  the  components  of  R°;  ©,  <X>  are 

its  coordinates  in  degrees-  Thus  X  =  R°  sin  9  cos  <t>,  Y  =  R°sln  ©  sin  *, 
o 

Z  =  R  cos  9 .  The  techniques  converge  very  rapidly  for  these  values  of 
K.  Even  for  Sample  3*  only  3  iterations  are  needed  to  obtain  <  =  3-92. 

4.  TESTS  OF  HYPOTHESES. 

4.1  We  now  consider  tests  for  bimodnl  data.  The  tests  to  be 
proposed  are  devised  to  make  use  of  methods  and  tables  already  prepared 
for  the  Fisher  distributions;  these  may  be  briefly  summarized  as  follows. 

For  theFisher  distribution,  one-sample  tests  of  hypotheses  concer¬ 
ning  A  or  K  are  based  on  R  and  on  X  (Stephens,  1962,  1967).  For 
the  important  test  that  s  different  samples  have  the  same  modal  vector 
A,  a  test  is  based  on  the  conditional  distribution  of  R^+I^+’-.+R  , 

given  R;  R  is  the  length  of  R,  the  overall  resultant  of  all  the  s 

samples.  The  tables  for  this  test  are  in  Stephens  (1969). 

4.2  Possible  test  statistics  for  the  Blmodal  Distribution:  R°  and  S. 

For  the  bimodal  distribution,  when  A  is  not  known,  it  would  be 

natural  to  base  tests  on  the  distribution  of  R°;  when  A  is  known, 
possible  test  statistics  could  be  the  component  of  R°  on  A,  or  the 
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sum  C  of  the  components  on  A  of  all  the  given  vectors,  each  one 
pointing  so  that  its  component  is  positive.  (Note  that  these  two 
statistics  are  not  necessarily  the  same.)  Unfortunately,  even  the 
distribution  of  R°  is  not  known.,  and  we  consider  an  alternative 
statistic.  Suppose,  for  a  given  sample  of  size  N,  we  have  found  the 

A  A 

estimates  A  and  <  by  the  above  methods;  we  then  ask  for  the  resul¬ 
tant  vector  g  which  would  have  given  the  same  estimates,  of  modal 
vector  A  and  of  JC,  on  the  assumpti-"  i  it  is  the  resultant  of  a 
Fisher-distributed  sample  of  the  same  size.  Clearly  g  will  lie  along 

R°,  and  its  length  is  easily  found  from  the  calculations  leading  to 

A 

K  for  the  bimodal  sample.  Consider,  for  example.  Sample  1.  The  final 
K  is  5-56,  obtained  by  solving  coth  <  -  1 Jk  =  0.8200:-  comparison  with 
equation  (2)  for  the  Fisher  distribution  shows  that  0.8200  must  he  the 
value  of  S/N,  so  that  S  =  0.8200  x  25  =  20-50.  We  call  S  the 
adjusted  resultant  of  the  bimodal  sample;  note  that  it  might  point  along 
either  direction  of  the  estimated  modal  axis,  according  to  the  final 
direction  taken  by  R°.  in  effect  we  have  imagined  constructing  a 
Fisher  sample  from  the  original  data,  with  same  N,  K}  and  modal  vector 
(based  on  our  estimates),  and  S  would  then  be  its  resultant.  The 
values  of  S  are  included  in  Table  2.  Since  we  hare  tests  and  tables 
available  for  the  Fisher  distribution,  based  on  the  Fisher  resultant 
R,  we  could  now,  as  an  approximate  procedure  for  the  bimodal  distribu¬ 
tion,  use  the  same  tests,  with  the  at ?d  resultant  g  replacing  g. 

One  sample  tests.  This  procedure  will  first  be  illustrated 


with  one-sample  tests  on  Sample  1. 


Test  for  the  modal  axis  A  Suppose  the  null  hypothesis  is  Hq: 

for  Sample  1,  is  along  $  =  90,  i*e*  the  y-axis.  The  appropriate 

Fisher  distribution  test  is  based  on  the  conditional  distribution  of 


R  given  the  component  C  or  the  hypothesized  modal  vector*  if  R  is 
too  large,  Hq  is  rejected.  Here,  the  value  of  S  is  20.50*  it 
lies  along  the  direction  of  R°,  and  its  component  C  is  then  18.22. 
The  test,  for  Q  =  0.05,  uses  Figure  2  of  Stephens  (1962),  or  approxima¬ 
tions  which  accompany  t  for  the  case  when  the  component,  called  X  in 


the  Figure,  is  beyond  the  range  given.  In  this  case,  the  approximation 


in  section  3.3  applies*  the  critical  value  of  S,  say  SQ,  will  be 


calculated  from 


So-°  F2,48(a> 


*  F2  48^°^  is  tte  usual  ^-distribution, 


here  with  2  ar.d  48  degrees  of  freedom,  at  the  level  ct.  For  Q  =  C.C5, 


we  get  S  =  19.016.  Since  S  is  greater  than  S  .  is  rejected 
O  O'  v 

at  the  5$  level. 


Use  of  R°.  For  reasonably  large  K,  R°  very  nearly  equals  g,  as 

here,  and  the  test  could  be  made  using  R°  }  the  component  C  is  now 

18.26  (the  Y  component  of  R  ,  from  Table  2)  and  the  critical  value 

for  R°  would  thor.  be  19-051.  H  would  again  be  rejected.  In  the 

o 

next  test,  for  K,  replacement  of  S  by  R°  al^o  does  not  change  the 


result  of  the  test. 


Test  for  x.  Consider  a  5$  test  of  Hq:  the  true  k  of  Sample  1  is 
4.5.  The  test  for  the  Fisher  distribution  is  given  in  Stephens  (1969), 
it  uses  the  value  of  I^/n,  and  a  table  of  the  critical  values  is  given- 
For  N  =  25,  x  -  4-5,  the  upper  5#  critical  value  is  0.854.  If  for 
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Sample  1,  we  use  R°  ae  tc  at  statistic,  the  value  of  R°/N  =  0.822; 
this  is  not  8lgnlf leapt,  so  is  not  rejected,  for  a  one-tail 

(or  two-tall)  test.  Similarly^  C,  N  =  0.820  and  use  of  S  would  give 
the  same  result. 


4.4  Tests  for  f-s/eral  samples.  Notation. 


Tut  the  r-oh  sample  have  size  N  ,  and  let  R  be  the  resultant 

r  ~r 

as  calculated  above,  and  S^,  the  adjus"  ed  resultant,  with  length  R°, 

S  i  as  before,  N  =  it  +  N.  •  •+  R  ,  where  s  is  the  number  of  samples, 

r  12, 

We  can  combine  the  resultants  P°  or  5  ,  or  the  samples  themselves, 

■  ~r 

i  ,  / 

in  several  ways: 

i  i 

(a)  Firstly,  fcupposu  vbe  rtsnl+ants  of  the  individual  samples 
j,  k,  1,  m,  say,  are  so  aligns:  that  they  give  the  maximum  length 

O  Q  O  O 

to  the  \ector  sum  R.  ■»  R,  -*  R  +  R  ;  we  shall  describe  them  as  well- 
aligned.  (Recall  that  any  cal  .ulated  R°  is  ambiguous  in  direction 
and  may  be  reversed  as  desired.)  For  well-aligned  vectors,  this  vector 

O  0 

sum  will  be  called  R.,  .  >  its  length  is  R..  , 

~jklm  a  jklm 

(b)  Secondly,  samples  j,  k,  1,  m,  may  all  be  pooled  into  one 
sample,  and  the  procedure  above  followed  to  calculate  R°  for  the 

o 

overall  sample)  we  shall  use  the  notation  P.,  ,  for  the  resultant 

~jklm 

of  pooled  samples.  It  will  not  necessarily  be  the  case  that 

will  be  equal  to  ;  a  vector  may  point  one  way  in  sample  j,  say, 

o 

in  calculating  R^,  and  be  reversed  in  calculating  the  pooled  sample 
resultant . 
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This  car.  easily  be  checked  by  reference  to  the  individual  components 

the  final  X  component  of  the  pooled  sample  will  equal  the  sum  of  the 

individual  X  components  if  the  vector  sum  equals  the  pooled  resultant, 

and  similarly  for  the  Y  and  Z  components.  As  an  illustration, 

consider  Table  2.  Samples  3  and  4;  even  whan  £°  is  reversed,  to  make 

tne  X  values  of  3  and  4  both  positive,  their  sum  is  not  the  X  value 

o 

of  Samples  3  and  4  pooled  together;  we  can  see  that  the  X  of  %4 
is  28-99  (16.74  +  12.25),  while  that  of  P°^  is  28. 16,  in  this  case 
there  occurs  one  vector  in  Sample  3  which  changes  direction  in  the 
overall  pooled  sample;  it  is  at  approximately  89.3  degrees  to 
ar.d  at  10?. 6  degrees  to  P°^.  l'n  such  circumstances  P°k  will  always 
be  greater  than  R°^  ,  has  Length  =  38.20  and  1°^  has 


length  P,jt<  =  38.62.  For  more  than  two  samples  the  disparity  may  be 

greater.  In  a  similar  way  w'  can  define  S,,  ,  =  S  +  S,  +  S,  ♦  S 

■'Jklm  ~k  ~1  ~m 

when  the  S  vectors  ar<=>  well-aligned;  its  length  is  S.,  ,  •  It  will 

jklm 

not  necessarily  be  the  same  as  the  value  of  S  for  the  pooled  samples 

k,  1,  and  mi  we  snail  call  this  ,  ,  length  Q.,  , 

‘  Kim  jklm 


4.c  Test  Statistics.  Tests  for  several  samples  might  be  based 

or.  R  or  S  ataM s1  ics.  and  or.  vector  sums  or  pooled  values-  We 

shall  suggest  using  S- statistics,  since  the  appropriate  tables  to  be 

used  are  based  or.  Fisher  distributions,  and  pooled  resultant  S 

values  (Q.,  .  rather  than  S.,  ,  ).  since  these  are  easily  obtained 
jslm  jAim 

from  the  computational  procedure  described  above,  when  the  samples  are 
pooled  into  ere.  Borderline  decisions  will  in  any  case  be  treated 


with  reserve  owing  to  the  approximate  nature  of  the  testa.  We  now 
illustrate  several  multi-sample  test  procedures. 

Test  for  a  Common  Modal  Axis.  A  test  for  Hq:  that  s  samples  with 
the  same  unknown  <,  have  a  common  modal  axis,  will  be  based  on  the 
arithmetic  sum  si+S2  +  *‘‘  +  ss  Siven  Q^g  .  s“  Suppose  we  test  thiB 
hypothesis,  at  the  5$  level  for  samples  1  and  2.  S1  =  20-50,  Sg  =  20-78, 
and  Qjg  =  4o.47«  The  test  follows  the  procedure  in  Stephens  (19&9)- 
The  steps  are  as  follows: 

(a)  Calculate  W  =  Q^N  =  0.8l0. 

(b)  Calculate  Z  =  (S^+SgVN  =  0-826. 

(c)  Use  section  2  of  Stephens  (1969)  to  find  the  critical  value 
z  for  N  =  50,  W  =  0.8l0,  and  Ot  =  0.05  i  using  the  F-approximation 
there  given,  we  find  z  -  0.821. 

(d)  Since  Z  exceeds  z,  reject  Hq  at  the  5$  level.  Z  is 

in  fact  Just  significant  at  the  2.5)6  level  (critical  value  z  =  0.824). 

The  vector  sum  S^g  =  Si+§2  bas  length  S^g  =  4o.49j  obviously 
replacement  of  Q^g  by  S-jg  gives  the  same  result  for  the  test. 

If  one  uses  the  R  statistics,  the  pooled  resultant  P^g  =  R-^g.  T^ie 
steps  (a)  and  (b)  give  W  =  0.811  and  Z  =  O.827.  Again  Z  is  just 
significant  at  the  2.5)6  levelj  the  critical  values  go  up  very  slightly 
with  the  slight  increase  in  W,  but  so  also  does  Z.  The  example 
illustrates  again  that  for  reasonably  large  K,  the  S  and  R  statis¬ 
tics  will  give  the  sane  results.  This  is  especially  so  in  the  tonts  of 
a  conditional  nature:  those  for  the  modal  axis  for  one  sample,  and  those 
for  a  common  modal  axis  for  several  samples . 
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5.  THE  VON  MISES  DISTRIBUTION  AND  BIMODAL  EXTENSION. 

5 . 1  The  distribution  analogous  to  Fisher'  b,  for  use  in  two  dimen¬ 
sions,  is  the  von  Mises  distribution;  if  9  is  the  polar  coordinate  of 
P,  or  equivalently  of  OP,  the  density  is 

f(6)  =  g"f  exp(*  cos  9)  -n  <  0  <  n 

with  modal  direction  A  along  9=0  and  k  measuring  concentration 
as  before.  Treatment  of  this  distribution  is  in  Gumbel,  Greenwood  and 
Durand  (1957)  and  Watson  and  Williams  (1956);  tests  are  in  Stephens 
(1962a,  1969b);  Batschelet  (1965)  gives  a  good  account  of  the  statistical 
procedures  and  a  bibliography  of  applications,  and  May  (1967)  shows 
how  the  procedures  may  be  applied  to  practical  situations.  Again 
estimation  and  testing  procedures  are  based  on  the  sample  resultant  R 
and  on  X,  its  component  along  A,  known  or  hypothesized.  The  direc¬ 
tion  of  g  estimates  A  when  this  is  not  known;  K  is  estimated  from 
Il(k)/1o(k)  =  R/N  ;  when  A  is  known,  X  replaces  R.  IQ(<),  Ip(*c) 
are  the  usual  imaginary  Bessel  functions  of  order  zero  and  one. 


5  .2  The  Blmodal  Extension.  The  extension  of  the  von  Mises  distri¬ 
bution  for  bimodal  data  gives  density 

(7)  f(9)  =  (a  exp(<  cos  9)+{  l-a)exp( -rx  cos  9))  ,  -#  <9_<_  «  ; 

o' 

when  a  =  0.5,  this  becomes 
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(8) 


t(9)  =  cosh(*  cos  0)  ,  -n  <  0  <  «  } 

with  &  modal  diameter  4  along  9=0  or  «.  If  &  lies  along 

9  »  9  or  0  +  it  ,  cos  9  is  replaced  by  cos(9-0  )  in  (7)  or  (8), 

0  0  o 

For  distribution  (8),  used  for  axial  data  or  directed  data  with 
opposite  inodes  of  equal  amplitude,  the  technique  for  estimating  A 
from  a  Bample  of  size  N  is  similar  to  that  in  section  3  for  three 
dimensions*  A  suitable  set  of  rectangular  coordinates  puts  the  x-axis 
along  0=0  (now  no  longer  known  to  be  the  modal  diameter ^A),  so 
that  x  =  cos  0,  y  =  sin  0  . 

5  *3  Estimation  of  R*  Then,  with  suitable  initial  unit  vector  v, 
components  1,  m,  A  is  estimated  iteratively  by  following  steps  (c) 
to  (e)  of  section  3*2  with  the  obvious  change  to  two  dimensions*  If 
no  good  estimate  of  A  exists,  we  follow  (f),  starting  now  with 
(l,m)  ■  (1,0)  and  then  (0,l). 

5*4  Estimation  of  < ■  The  equation  for  estimating  k  is  now 

l  N 

(9)  I1(/c)/lo(x)  =  X  cos  0jtanh(K  cos  0^ 

and  this  is  solved  iteratively  exactly  as  in  section  3*3,  steps  (a) 
to  (c)>  the  initial  right  hand  side  is  R°/N,  where  R°  is  the  length 
of  the  vector  estimating  A.  The  successive  k  values  must  be  found 
by  interpolation  in  a  table  giving  values  of  11(k)/Iq{k)  for  given 
a  table  is  in  Gumbel,  Greenwood  and  Durand  (1953)  or  Batschelet  (1965)* 
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5 .5  Examples .  Some  axial  data,  supplied  by  Dr.  A.  Rees  of  the  Depart 
ment  of  Oceanography,  University  of  Southhampton,  are  in  Table  3, 

Samples  5  and  6.  The  data  represents  axis  of  maximum  susceptibility  in 
magnetization  of  rocks,  and  comes  from  the  Franciscan  rockB  in  Diablo 
Range  in  California,  (Rees  and  Hamilton,  1965)-  Table  3  gives  the 
estimates  of  the  parameters  (S,  the  polar  coordinate  of  JJ,  gives 
the  estimated  inclination  of  A),  for  the  samples  taken  separately 
and  also  pooled  into  one  sample.  The  modified  statistic  jg  is  calculated 
as  for  three  dimensions ;  the  final  right-hand  side  of  (9)  gives  S/M, 
and  S  lies  along  R°.  Tests  are  conducted  as  for  three  dimensions; 
those  for  the  modal  vector  and  for  <  are  in  Stephens  (1962a,  1962b); 
multi-sample  tests  are  in  Stephens  (1969c).  We  willustrate  only  a 
two- sample  test  with  Samples  5  and  6. 

Test  for  a  Common  Modal  Axis.  For  a  test  of  the  two  samples  have 

a  common  modal  axis  we  have  Z  =  (S^  +  Sg)/N  -  0-973  ;  W  =  Q^^/N  =  0.959* 

The  critical  value  z  of  Z  is  given  in  Stephens  ( 1969c);  this  is  a 
slight  improvement,  or.  an  F-test  introduced  by  Watson  and  Williams 
(  1956).,  and  gives  z  =  C.986  at  the  5$  level-  Thus  we  reject  the  null 
hypothesis  Hq  at.  this  level. 

6.  OFPOSITE  MODES  OF  UNEQUAL  AMPLITUDE- 

6.1  Data.,  with  direction,  from  populations  with  modes  which  are 
opposite  but  of  unequal  amplitude,  will  be  treated  using  distributions 
(3)  or  (7)  The  vectors  now  have  a  definite  sense.  The  estimation 
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of  A  follows  the  steps  as  before,  in  sections  3  and  5}  then  a  may  be 
estimated  by  a  ,  the  proportion  of  vectors  which  is  not  reversed  when 
estimating  A  .  The  equations  for  estimating  K  become 

1  1  N  a  exjj(K  cos61)-(l-a)exp(-K  cosS^ 

(10)  coth  ft  --  =  5  Z^oa  01  —  (K  c08ei)+(1-a)exp(-K  coTeJ 

for  three  dimensions}  for  two  dimensions  the  left  hand  side  of  (10)  is 
replaced  by  I-L(x  )/lo(/c)  •  In  the  right  hand  side  one  inserts  a  ,  and 
follows  the  same  iteration  procedure  as  in  section  3}  the  sequence  for 
K  converges  to  the  estimate  k  • 

Example •  We  illustrate  this  technique  with  some  interesting  data  given 
several  years  ago  by  Dr.  E.  Gould  of  the  Johns  Hopkins  School  of  Hygiene. 
The  data  represent  directions  taken  by  turtles  after  treatment;  it  is 
thought  that  the  turtles  have  a  preferred  direction,  but  some  are 
confusing  forwards  with  backwards.  Thus  the  distribution  is  (7)}  the 
actual  values  and  the  analysis  are  in  Table  3>  May  (1967)  has  varied 
the  parameters  to  attempt  to  find  a  best  fit;  his  best  fit  values  are 
9  =  61.5,  K  =  3.167,  and  a  =  0.803 •  The  results  in  Table  3  are  in 
excellent  agreement . 

Tests  and  Confidence  Intervals-  These  would  follow  the  procedures 
already  described,  using  again  the  right  hand  side  of  (10)  to  give  an 
adjusted  resultant  S.  In  this  example  it  would  have  length  61.83,  and, 
ub  before,  lies  along  0  =  63.08,  the  direction  of  g.  Suppose  one  wished 
to  find  a  1056  confidence  interval  for  the  modal  axis  of  (7),  using  g; 


Figure  1  of  Stephens  (1962a)  is  used,  or  the  approximations  given  when 
the  data  is  beyond  the  range  of  the  figure}  here,  the  approximation  in 
section  3 .4  applies,  and  gives  a  band  for  9  with  a  half -width  equal 
to  7.4  degrees.  The  confidence  interval  is  55-7  <  ®  <  70*5  • 

7.  FURTHER  REMARKS. 

7.1  importance  of  knowledge  of  type  of  data.  We  illustrate 
this  point,  mentioned  in  section  3-1,  with  another  two-dimensional 
sample  of  Rees  and  Hamilton  (1965);  the  data,  io  values  of  9,  1b 
from  their  site  San  Jose  9-  Stoe  values  are  2,  13,  14,  l4l,  152,  156,  166 

356,357,358-  A  rough  glance  might  suggest  a  von  Mises  sample;  if  so 
analyzed,  the  modal  direction  would  be  along  0  =  41.13;  the  resultant 
is  3- 18;  K  is  0.672,  indicating,  of  course,  vide  dispersion.  In 
fact  the  data  is  axial;  if  diameters  are  drawn  through  the  data  points, 
we  see  one  end  of  each  diameter  makes  a  set  concentrated  between  260 
and  20;  if  the  analysis  of  this  paper  is  used,  the  modal  axis  is 
3  =  351*6,  R°  is  S-59,  S  =  9*59,  K  =  12.10.  There  is  a  big 
difference  in  the  estimates  of  modal  direction.  A  knowledgeable  ex¬ 
perimenter  may  of  course  present  the  data  so  that  either  analysis  would 
give  the  same  result  for  modal  direction;  this  has  happened  in  our  data 
in  Sample  1,  where  at  the  end  we  see  that  R°  is  found  with  none  of  the 
original  vectors,  as  given,  needing  to  be  reversed.  Thus  straight 
Fisher-distribution  techniques  would  have  given  the  same  resultant. 

But  it  would  have  been  easy  to  reverse  a  selection  of  the  given  vectors 
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to  produce  an  entirely  different  Fisher  resultant.  On  the  sphere, 
especially,  it  is  not  easy  to  see,  especially  with  widespread  data, 
which  end  should  be  chosen  so  that  a  Fisher  analysis  gives  the  same 
resultant  as  the  techniques  of  section  2}  the  point  of  these  techniques 
is  to  render  such  a  choice  unnecessary.  Note  that  in  any  case,  even 
with  the  same  resultant,  one  would  obtain  a  different  <  for  Fisher- 
distribution  and  for  bi-modal  analysis. 

7.2  Alternate  distributions-  Another  distribution  for  bimodal 
or  axial  data  has  been  proposed  by  Watson  ( 1966 ) >  the  density  for  0 
is  proportional  to  exp(X  coa^0)  (if  the  modal  axis  A  is  along 
0  «s  0),  and  ®  is  uniform  as  before.  Use  of  this  distribution,  with 
estimation  as  described  in  Watson  (1966)  gives  the  following  modal  axis  estimates 
for  Samples  1  to  4  in  Table  1)  Sample  1,  0  =  102.6,  <t>  =  115 .6  » 

Sample  2,  0  =  78-2,  •  =  111-7  i  Sample  5,  0  =  98-2,  *  =  148.4  >  Sample 
4,  0  =  106.2,  ®  ■  127.9.  The  results  are  in  good  agreement  with  those 
given  in  Table  1.  Testing  procedures  are  not  yet  as  well  developed 
for  this  distribution  as  for  the  Fisher  distribution. 

In  two  dimensions,  the  corresponding  density  is  equivalent  to  a 
density  proportional  to  exp(l  cos  20 )j  thus  the  doubled  angles  have  a 
von  Mises  distribution,  and  analysis  proceeds  by  doubling  the  angles 
given,  estimating  the  modal  vector,  and  halving  its  angle-  Thus  for 
the  doubled-angle  vectors,  for  Sample  5,  R  =  7 *40,  along  0  =  185 .6, 
so  that  the  modal  axis  estimate  is  A  ,  along  0  =  91*8  »  for  Sample 

A 

6,  R  is  8.76,  along  6  =  222.5,  giving  modal  axi6  A  along  0  =  111.25- 


22 


Direct  application  of  voc  Mises  techniques  rejects  the  hypothesis.,  at 
the  5#  level,  that  the  samples  have  a  common  modal  vectorj  details  are 
in  Stephens  (1969c).  These  results  compare  well  with  those  in  Table  3* 
Neither  of  the  above  distributions  in  two  or  three  dimensions,  is 
strictly  applicable  to  directed  data  with  unequal  modes,  except  that  tbc 
same  techniques,  to  estimate  the  modal  axis,  cams  till  be  employed.  If 
this  is  done  with  the  turtle  data  of  Sample  7,  i.e.  the  angles  are 
doubled  and  the  direction  of  the  resultant  then  found  ana  halved,  we 
have  the  estimate  of  A  along  0  =  62.57- 

7-3  Goodness-of -f It .  One  might  wish  to  test  the  data  to  see  if 
they  are  well-fitted  by  the  distributions  considered.  Trie  is  an 
important  subject,  and  for  the  present,  only  a  few  comments  are  offered 
here.  For  both  sphere  and  circle,  a  distinction  must  be  made  between 
axial  and  timodal  directed  data.  When  the  axis  A  has  been  estimated, 
it  should  be  taken  as  origin  for  9  and  all  coordinates  transformed. 
Axial  data  should  be  so  recorded  that  all  vectors  are  within  90  degrees 

^  JT 

of  A,  i  e.  their  0-values  are  less  than  rr  radians-.  For  the  sphere. 

the  test  is  then  made  for  a  fit  to  distribution  (5)-  .  For  bimodal 

directed  data,  (4)  is  used,  or,  with  unequal  inodes,  (3)  The  tests  can 

2 

be  made  separately  for  9  and  for  <5>,  using  say  the  X  test.  A 
rough  measure  of  goodness  of  fit  can  be  found  by  -use  of  the  U2  and  V 
goodness-of -fit  statistics,  but  no  precise  tests  can  be  made/  as  their 
distributions  depend  or.  the  fact  that  parameters  have  been  estimated. 
Similar  remarks  apply  to  the  bimodal  distributions  on  the  circle.  When 


the  coordinates  he.ve  been  transformed  to  make  6=0  lie  along  A, 
distributions  (7)  or  (8)  are  used  for  bimodal  directed  data.  For  axial 
data,  distribution  (8)  is  used,  with  twice  the  f(0)  shown  and  range 
-«/2  <  6  <  jc/2  . 
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TABLE  1 


Three-dimensional  axial  data:  estimation  procedures  for  Sample  1 
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Sample  1:  Coordinates  0,  $  of  25  vectors:  1  Is  the  vector  number. 
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Estimation  of  modal  axis  A»  final  estimate-  underlined. 
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.rj-;  final  resultant. 


20.54;  adjusted  Fisher  resultant  S.  -  20.50 
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